-- https://rosettacode.org/wiki/Euler_method#Common_Lisp
--
-- Specialised to the cooling function.  We produce an array of the
-- temperature at each step subtracted from the analytically
-- determined temperature (so we are computing the error).
--
-- ==
-- input { 100.0 0.0 100.0 2.0 }
-- output {
--   [0.0f64, 0.7486588319044785f64, 1.2946993164580505f64, 1.679265585204547f64,
--    1.9360723079051994f64, 2.0926628953127633f64, 2.1714630634463745f64,
--    2.1906621307551717f64, 2.1649494013163135f64, 2.106128735048543f64,
--    2.0236308042860642f64, 1.924939486661728f64, 1.8159462635058006f64,
--    1.7012443098140402f64, 1.5843721158726787f64, 1.4680149205211634f64,
--    1.354170918183268f64, 1.2442880891724784f64, 1.139376563939308f64,
--    1.0401006401924064f64, 0.9468539045506148f64, 0.859820348366565f64,
--    0.7790238943202752f64, 0.7043683525135087f64, 0.6356694904077109f64,
--    0.572680620138101f64, 0.5151128711182906f64, 0.46265111832983763f64,
--    0.4149663712668854f64, 0.37172529010454625f64, 0.3325973799923041f64,
--    0.29726031781953566f64, 0.2654037853016753f64, 0.23673211521330373f64,
--    0.21096600187360082f64, 0.1878434807373779f64, 0.16712034361770378f64,
--    0.14857012436523576f64, 0.13198376366226583f64, 0.11716904003828077f64,
--    0.10394983650827427f64, 9.216529772764304e-2f64, 8.166892070760667e-2f64,
--    7.232761248797459e-2f64, 6.402074034058458e-2f64, 5.663919376000237e-2f64,
--    5.008447242299141e-2f64, 4.426781024452353e-2f64, 3.9109342442220196e-2f64,
--    3.453731999022125e-2f64]
--
-- }
--
-- input { 100.0 0.0 100.0 5.0 }
-- output {
-- [0.0f64, 4.375047177497066f64, 5.926824303312763f64, 6.025019928892426f64,
-- 5.447257115328519f64, 4.619590476035608f64, 3.763003010238556f64,
-- 2.9817046074496396f64, 2.315646506892435f64, 1.7712171223319615f64,
-- 1.338771206215167f64, 1.0023162611494705f64, 0.7446054205897603f64,
-- 0.5495998794267152f64, 0.4034719393409745f64, 0.2948359128180762f64,
-- 0.21460148389295242f64, 0.15566929274598706f64, 0.11258571548996343f64,
-- 8.121463510024185e-2f64]
-- }
--
-- input { 100.0 0.0 100.0 10.0 }
-- output {
-- [0.0f64, 15.726824303312767f64, 12.52775711532852f64, 7.636514260238556f64,
-- 4.216805010017435f64, 2.221390673785482f64, 1.1413261456382173f64,
-- 0.5782306456739441f64, 0.290580297318634f64, 0.14532974216231054f64]
-- }

def analytic (t0: f64) (time: f64) : f64 =
  20.0 + (t0 - 20.0) * f64.exp (-0.07 * time)

def cooling (_time: f64) (temperature: f64) : f64 =
  -0.07 * (temperature - 20.0)

def main (t0: f64) (a: f64) (b: f64) (h: f64) : []f64 =
  let steps = i64.f64 ((b - a) / h)
  let temps = replicate steps 0.0
  let (_, temps) =
    loop (t, temps) = (t0, temps)
    for i < steps do
      let x = a + f64.i64 i * h
      let temps[i] = f64.abs (t - analytic t0 x)
      in ( t + h * cooling x t
         , temps
         )
  in temps
